DECOMPOSITION  OF  SUPERIMPOSED  WAVEFORMS  USING  THE  CROSS 

TIME  FREQUENCY  TRANSFORM 


P.  Bonato,  Z.  Erim,  J.A.  Gonzalez-Cueto 
NeuroMuscular  Research  Center,  Boston  University,  Boston,  MA  USA 


Abstract  —  The  identification  of  the  timing  of  the 
discharges  of  groups  of  muscle  fibers  (motor  units)  is  of 
utmost  importance  in  research  into  the  strategies  employed  hy 
the  central  nervous  system  in  producing  muscle  force  and  in 
the  diagnosis  of  neuromuscular  diseases.  The  process  involves 
the  recognition  of  unique  shapes  (action  potentials) 
contributed  by  different  motor  units  at  random  times 
throughout  a  muscle  contraction.  This  paper  addresses  a 
specific  aspect  of  the  identification  process:  the  decomposition 
of  the  compound  signal  when  the  action  potentials  of  two  or 
more  motor  units  are  superimposed.  We  propose  a  cross-time- 
frequency-based  procedure  to  identify  which  two  (out  of  a 
previously  identified  collection  of  waveforms)  are  included  in 
a  superposition. 
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I.  INTRODUCTION 


waveforms,  referred  to  as  ‘templates’,  have  already  been 
identified.  Therefore  let  Xj ,  X2 ,  ...  ,  be  the  templates 
derived  for  N  motor  units  and  y  a  signal  constituted  by  the 
superposition  of  and  X j  ,which  needs  to  be  decomposed.  The 

technique  that  we  propose  consists  of  deriving  the  cross-time- 
frequency  distribution  [6]  of  each  template  with  the  signal  y 
constituted  by  the  superposition  of  any  combination  of  the 
identified  action  potential  waveforms  with  unknown  delays  with 
respect  to  each  other.  It  can  be  shown  that  when  the  template 
utilized  to  derive  the  cross-time-frequency  distribution  is 
embedded  in  the  signal  y,  the  cross-time-frequency  distribution  - 
filtered  according  to  criteria  designed  in  the  ambiguity 
domain  [10]  -  is  “similar”  to  the  time-frequency  distribution  of  the 
template.  On  the  contrary,  when  the  template  utilized  to  derive  the 
cross-time- frequency  distribution  is  not  embedded  in  the  signal  y, 
the  filtered  cross-time-frequency  distribution  is  not  “similar”  to 
the  time- frequency  distribution  of  the  template. 


A  variety  of  algorithms  have  been  proposed  in  the  past  to 
decompose  the  signal  into  the  basic  waveforms  related  to  the 
firing  of  different  motor  units.  However,  when  two  or  more 
different  action  potential  waveforms  occur  at  the  same  time,  the 
automatic  decomposition  of  the  signal  into  the  two  action 
potential  waveforms  is  not  a  trivial  task.  Traditional  filtering 
approaches  are  not  very  successful  in  this  case  because  motor  unit 
action  potential  waveforms  often  partially  overlap  in  frequency. 
Therefore  trivial  approaches  (i.e.  matched  filters)  are 
unsatisfactory  to  solve  this  problem  [8].  The  most  common 
approach  previously  proposed  is  the  peel-off  method  [4][11] 
where  the  superposition  is  correlated  to  each  individual  motor  unit 
template  and  the  best  matched  template  is  aligned  with  the 
superposition  and  subtracted  from  it.  The  procedure  is  repeated  for 
the  rest  of  the  templates  until  the  energy  of  the  remaining  signal  is 
below  a  preset  value.  This  approach  is  also  unsatisfactory  because 
depending  on  the  delay  with  which  the  action  potentials  overlay 
each  other  various  peaks  can  cancel  each  other  out  and  the 
resultant  signal  may  not  have  strong  correlation  to  the  individual 
templates  [8].  In  order  to  overcome  the  limitations  of  previous 
approaches  we  propose  to  investigate  a  methodology  that  relies  on 
a  cross-time-frequency  based  algorithm.  The  cross-time-frequency 
representations  will  be  computed  using  Cohen  Class 
transformations  [7]  because  of  the  properties  of  this  class  of 
transformations  that  allow  the  filtering  out  of  the  cross-terms 
[9][12]  that  make  it  difficult  to  solve  the  superposition  problem. 

II.  CROSS  TIME  FREQUENCY  BASED 
DECOMPOSITION  OF  SUPERPOSITIONS 

A.  Rationale 

The  proposed  methodology  to  decompose  superpositions 
will  assume  that  a  certain  number  of  unique  action  potential 


B.  Illustration  of  the  Technique 


In  the  following  we  illustrate  the  procedure  for  a 
superposition  =  Xj  -F  Xj  .  We  first  consider  the  cross-time- 
frequency  representation  of  Xj  and  y ,  namely  ,  and 

observe  that  it  can  be  expressed  as 


TF  =TF  +  TF 

XiV  x,Xi  X, 


where  TF^^^^  is  the  time-frequency  representation  of  Xj , 
and  77c  is  “half’  the  interference  term  [2]  of  the  time- 

X1X2  L  j 

frequency  representation  of  y ,  TF^^  expressed  by 
TFyy  =  TF^^^^  +TF^^^^  +^^X,X2  +^‘^X2X1-  By  applying  an 
appropriate  kernel  to  the  cross-time-frequency  transform  TF^  ^  , 
it  is  possible  to  reject  TF^^^^  thus  leading  to  TF^^^^  . 

This  is  because  77c  ,  holds  the  characteristics  of  the 

X1X2 

interference  terms  [3]  and  thus,  in  order  to  reject  it,  one  may  apply 
the  same  criteria  utilized  to  attenuate  the  interference  terms  which 
affect  the  auto-time- frequency  representation  of  a  generic  signal. 
In  practice  an  attenuation  of  77c  ,  will  be  obtained  rather  than 

X  X,X2 

its  complete  rejection  [12]. 

The  technique  is  illustrated  in  Figure  1.  The  upper  plot 
represents  the  action  potential  waveforms  of  two  templates,  Xi(t) 
and  X2(t).  The  plot  in  the  middle  shows  the  time- frequency 
representation  of  the  two  templates.  The  lower  plot  displays  the 
cross-time-frequency  representation  of  Xi(t)  and  y(t)=Xi(t)+X2(t) 
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Figure  1.  The  cross-time-frequency  based 
technique  to  solve  the  superposition  of  waveforms. 

(A)  Superimposed  waveforms  Xi(t)  and  X2(t). 

(B)  Time-frequency  distributions  and  . 

(C)  The  Choi-Williams  cross-time-frequency 
distribution  computed  using  the  first  template  Xi(t) 
with  the  superimposed  signal  y(t)=Xi(t)-l-X2(t). 


obtained  using  the  Choi-Williams  transformation  [5],  The  cross- 
time-frequency  representation  in  Figure  1  is  apparently  “similar” 
to  the  time-frequency  representation  of  Xi(t). 


+  AF, 


where  AF  indicates  either  the  ambiguity  function  or  the 
cross-ambiguity  function.  Subscripts  correspond  to  those  defined 
above  for  the  time-frequency  representations.  It  is  well  known  that 

the  term  AF  is  usually  concentrated  over  the  axes  of  the 
ambiguity  plane,  while  the  term  AF^  ^  ,  which  corresponds  to  a 
“half’  interference  term,  is  expected  to  be  away  from  the  axes. 
The  distance  of  AF^  ^  from  the  origin  of  the  ambiguity  plane  is 
a  function  of  the  distance  between  the  baricenters  of  the  auto- 
time-frequency  representations  of  Xj  ( TF^^^^  )  and  X2 

(  TF^  ^  ).  Therefore  kernels  which  attenuate  components  of  the 
ambiguity  function  that  are  away  from  the  axes  lead  to  a  reduction 
of  the  amplitude  of  TF^^^^  . 


Figure  2  shows  the  ambiguity  functions  of  the  action 
potential  waveforms  represented  in  Figure  1.  The  ambiguity 
function  of  Xi(t)  and  the  cross-ambiguity  function  of  Xi(t)  and  X2(t) 
are  displayed.  The  sum  of  these  two  functions  is  equal  to  the 
cross-ambiguity  function  of  Xi(t)  and  y(t)=Xi(t)  +  X2(t),  as  seen 
from  equation  (1)  above.  The  coordinates  of  the  baricenter  of  the 
cross-ambiguity  function  of  Xi(t)  and  X2(t)  are  equal  to  the 
distance  in  time  and  frequency  of  the  time-frequency 
representations  of  the  two  components  [9].  This  figure  also 
emphasizes  the  relevance  of  an  appropriate  choice  of  the  filtering 
function  in  the  ambiguity  domain.  The  optimal  filter  would 
maximize  the  attenuation  of  the  cross-ambiguity  function  of  Xi(t), 
X2(t)  and  not  alter  the  ambiguity  function  of  Xi(t).  Possible  design 
criteria  are  discussed  in  the  last  section  of  this  paper. 


D.  Classification  Procedure 


With  _y  —  X,-  -h  X j  and  i  and  j  being  unknown,  the 
proposed  classification  procedure  requires  the  computation  of  all 
the  cross-time-frequency  representations  of  X^  (with  k=l,  ...  ,  N) 
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C.  Interpretation  in  the  Ambiguity  Domain 

To  further  demonstrate  how  the  cross  time-frequency 
approach  suppresses  undesirable  terms  and  ‘teases  out’  the 
template  embedded  in  the  superposition,  consider  the  cross¬ 
ambiguity  function  derived  by  Xj  and  y 


Figure  2.  The  ambiguity  function  of  the  template 
X|(t)  and  cross-ambiguity  function  of  Xi(t)  and 
X2(t).  The  relationship  between  the  ambiguity  and 
cross-ambiguity  functions  are  used  to  filter  the 
cross-ambiguity  function  in  order  to  identify  the 
template  contributing  to  the  superposition. 


and  y  thus 


TF^  =  TF^ 

^ky 


+  TF^ 


and  the  cross-ambiguity  function 


AF  ^  =  AF  ^  ^  +  AF  ^  ^ 

^  k  y  ^  I  ^  k^  2 


It  is  apparent  that  in  cases  where  7>  i  or  j,  the  cross¬ 
ambiguity  function  is  expected  to  be  located  away  from  the  axes 
while  the  energy  around  the  origin  of  the  ambiguity  plane  is 
expected  to  be  negligible.  In  fact,  in  these  cases,  we  are  computing 
the  cross-ambiguity  function  of  cross-products  among  different 
components,  i.e.  templates,  thus  leading  to  terms  that  are  like  the 
interference  terms. 

An  example  is  shown  in  Figure  3  that  shows  the  output  of 
the  proposed  procedure  when  one  computes  the  cross-time- 


time  (ms) 


Figure  3.  The  rejection  of  a  template  that  is  not 
included  in  a  superposition.  (A)  Two  templates  Xi(t) 
and  X2(t),  that  contribute  to  a  superposition  y(t)=  Xi(t) 
+  X2(t),  and  a  third  template  X3(t)  that  is  not  included 
in  the  superposition.  (B)  Time-frequency 
representation  of  the  third  template  X3(t).  (C)  The 
Choi-Williams  cross-time-ffequency  distribution  of 
y(t)=  xi(t)  -f  X2(t)  and  X3(t). 


frequency  representation  of  the  superimposed  waveform 
y(t)=Xi(t)+X2(t)  with  a  template  Xjft).  The  upper  plot  of  Figure  3 
shows  the  motor  unit  action  potential  waveforms  Xi(t)  and  X2(t), 
that  contribute  to  the  superimposed  waveform,  together  with  a 
third  template  The  plot  in  the  middle  represents  the  time- 
frequency  distribution  of  the  template  X}(t).  The  lower  plot  reports 
the  cross-time-frequency  distribution  of  Xjffj  anA  y(t)=Xi(t)+X2(t) 
computed  using  a  Choi-Williams  transformation  [5].  Clearly,  the 
derived  cross-time-frequency  distribution  does  not  even  slightly 
resemble  the  time- frequency  distribution  of  Xi(t). 

Figures  1  and  3  summarize  the  rationale  for  the 
identification  procedure.  When  one  computes  the  cross-time- 
frequency  distribution  of  the  superimposed  waveform  and  a 
template  that  does  contribute  to  the  superposition,  the  computed 
representation  on  the  time-frequency  domain  is  “similar”  to  the 
distribution  of  the  template.  Conversely,  when  one  derives  the 
cross-time-frequency  distribution  of  the  superimposed  waveform 
and  a  template  that  does  not  contribute  to  the  superposition,  the 
derived  representation  on  the  time-frequency  domain  does  not 
even  resemble  the  time-frequency  distribution  of  the  template. 
Therefore,  using  a  measure  of  “similarity”,  i.e.  distance,  between 
the  cross-time-frequency  representation  and  the  time-frequency 
distribution  of  the  template,  one  may  define  a  way  to  identify  the 
templates  that  contribute  to  the  superimposed  waveform. 

E.  Filtering  of  the  Cross  Ambiguity  Function 

To  increase  the  “similarity”  between  time-frequency  and 
cross-time- frequency  representations  (Figure  2)  we  compared 
different  approaches  to  filtering  the  cross-ambiguity  function  (i.e., 
choices  of  the  kernel  of  the  cross-time-frequency  transformation). 
Specifically,  we  utilized  the  following  transformations:  1)  cross- 
Wigner-Ville  (XWV);  2)  Choi-Williams  (CW)  (o=0.1)  [10]; 
3)  Radially  Gaussian  Kernel  (RGK)  [1];  4)  a  transformation  based 
on  a  modified  design  of  the  RGK  kernel  (Modified  RGK),  derived 
by  utilizing  the  entire  ambiguity  function  domain  for  the 
convergence  of  the  algorithm  as  opposed  to  the  first  and  second 
quadrants;  5)  a  Matched  Kernel  transformation  designed,  for  each 
template  under  consideration,  as  the  square  of  the  magnitude  of  its 
ambiguity  function.  Figure  4  displays  the  improvement 
introduced  by  the  use  of  the  Matched  Kernel. 

F.  Simulations 

Five  actual  action  potential  waveforms  (templates)  recorded 
during  a  muscle  contraction  were  used  to  simulate  superimposed 
waveforms.  A  total  of  400  superpositions  were  obtained  by 
combining  pairs  of  templates  with  different  degrees  of  overlap, 
i.e.,  shifting  one  template  with  respect  to  the  other  (range  ±1  ms). 
For  each  superposition,  the  cross-time- frequency  representations 
of  the  superposition  and  each  of  the  templates  were  computed 
using  the  methods  described  above.  Then  the  Kolmogorov 
distance  between  each  cross-time-frequency  representation  and 
the  time-frequency  distribution  of  the  template  used  to  derive  the 
cross-time-frequency  representation  was  computed.  The 
superposition  was  considered  to  be  made  up  of  the  two  templates 
that  resulted  in  the  lowest  two  distances. 
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Figure  4.  Improvements  afforded  by  filtering  in  the 
ambiguity  domain.  A.  The  cross-ambiguity  functions 
between  the  waveforms  Xi  and  X3  with  the 
superposition  waveform  y  =  Xj  -I-  X2  (see  Figure  2). 
The  kernels  are  shown  as  shaded  areas.  B.  The 
resultant  cross-time-ffequency  transforms  between 
the  superposition  waveform  and  the  templates  Xj  and 
X3.  Note  the  increased  similarity  between  the  cross- 
time-ffequency  representation  and  the  time- 
frequency  distribution  of  Xj  included  in  the 
superposition. 

III.  PRELIMINARY  RESULTS 

Our  study  indicates  that  the  proposed  method  can  identify 
the  motor  unit  action  potential  waveforms  contributing  to  a 
superposition  with  reasonable  success.  Marked  differences  were 
observed  when  different  filtering  approaches  were  applied  to  the 
cross-ambiguity  functions.  The  XWV  approach  resulted  in  the 
lowest  performance,  i.e.,  58.5  %  of  the  decompositions  were 
successful.  The  Matched  Kernel  approach  was  the  best  (88  % 
successful  decompositions). 

IV.  CONCLUSIONS 

A  method  to  solve  the  superposition  of  action  potential 
waveforms  has  been  proposed  based  on  the  use  of  cross-time- 
frequency  transformations.  A  preliminary  study  to  assess  its 
application  to  real  data  indicates  the  suitability  of  the  technique. 
However,  future  work  needs  to  be  done  in  order  to  fully 
characterize  the  methodology.  Specifically,  a  point  of  paramount 
importance  is  the  characterization  of  the  methodology  when  the 
motor  unit  action  potential  waveforms  slightly  change  in  time,  as 
it  happens  because  of  the  physiology  of  the  muscle  contraction. 
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Figure  5.  Successful  decompositions  of  the  proposed 
procedure.  The  different  bars  show  the  results  for  the 
cross-Wigner-Ville  transform  (XWV),  the  cross- 
Choi-Williams  cross-transformation  (CW),  the 
radially  Gaussian  kernel  cross-distribution  (RGK),  its 
modified  version  (Modified  RGK),  and  the  Matched 
Kernel  approach. 
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